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I. 



INTRODUCTION 



A. DIGITAL SIGNAL PROCESSING 

Digital signal processing has its roots in 16th century 
mathematics, especially in the fields of astronomy and the 
compilation of mathematical tables. Today it has become 
a powerful tool in a multitude of diverse fields of science 
and technology. The applications of digital signal pro- 
cessing varies from low-frequency spectrum seismology 
through spectral analysis of speech and sonar into the 
video spectrum of radar systems [1] . 

In their book. Theory and Applications of Digital 
Signal Processing [2] , Rabiner and Gold have combined digi- 
tal signal processing theory, with a variety of applications 
ranging from sonar, radar, communication, music, seismic 
and medical signal processing, with digital component 
technology. This technology is the main driving force for 
progress in this field as well as the general area of com- 
puter design. It was also observed by them, that, although 
the formulation of engineering problems is often as vague 
as those of the "softer" sciences (such as anthropology, 
psychology, etc.), the execution of these problems appears 
to depend on greater and greater accuracy and reproducibility. 
The capability of digital systems to achieve a guaranteed 
accuracy and essentially perfect reproducibility is very 
appealing to engineers. 
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A digital filter is defined in [3] to be a computa- 
tional process in which a sequence of numbers acting as 
input, is transformed into a second sequence of numbers, 
or output digital signal; where the term "digital” implies 
that both time (the independent variable) and amplitude are 
quantized. 

The field of one-dimensional digital filtering, which 
is outlined in figure 1.1, encompasses recursive, non- 
recursive and Fast-Fourier Transform processing. The terms 
recursive and nonrecursive, instead of HR (Infinite Impulse 
Response) and FIR (Finite Impulse Response) are preferred 
in this thesis. It is shown that recursive processing is 
much more efficient than nonrecursive processing. Stockham's 
method [4] to perform fast convolution, which later became 
known as the FFT method, improved the efficiency of non- 
recursive techniques, so that comparisons in one^-dimension 
are no longer strongly biased toward recursive methods [ 2 ] . 

One-dimensional recursive digital filter design depends 
strongly on the effective and well developed continuous 
filter design theory. Moreover, the stability analysis of 
higher order one-dimensional recursive realizations can be 
solved by investigation of the root distribution of its 
factored form representation. It is important to realize 
that continuous domain design techniques and factorization 
property (the fundamental theorem of algebra) exists only 
in one-dimension. 
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FIGURE 1.1. OVERVIEW OF ONE-DIMENSIONAL DIGITAL SYSTEMS 




Areas of application of one-dimensional digital signal 
processing are listed in figure 1.1. There are important 
military and specifically naval applications in addition 
to those listed, such as missile and torpedo guidance, 
launcher control system, fire control, combat information 
collection, dissemination and display. 

Digital filtering in several dimensions, which is over- 
viewed in figure 1.2, has gained considerable importance, 
especially for the two-dimensional case. Figure 1.3 presents 
an example of two-dimensional low-pass and high-pass filtering. 
Until 1966, two-dimensional digital filtering was implemented 
by nonrecursive or convolutional techniques. In this method 
the output is the weighted sum of unit sample responses of 
all past input values. The serious disadvantage of the 
convolutional method is the requirements of a very large 
number of arithmetic operations. 

The development of the Fast Fourier Transform (FFT) in 
1966, reduced the number of arithmetic operations consider- 
ably and is used extensively today. Filtering via FFT is 
accomplished by computing the transform of the input func- 
tion, multiplying by the frequency response of the filter, 
and inverse transforming the result. The recursive 
algorithm, which has in general an infinite impulse response 
[3] , constitutes another technique for the realization of 
two-dimensional digital filters. 

Hall [5] compares the amount of computation required for 
the three filtering techniques. He shows that, in general, 
the FFT and recursive algorithms are preferable to the 
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FIGURE 1.2. OVERVIEW OF N- DIMENSIONAL DIGITAL SYSTEMS 



I 



Fig. 1.3 




.a. An example of two-dimensional filtering 
(An original photograph) [2] 
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Fig. 1.3.b. An example of two-dimensional filtering 
(Low pass filtered version of Fig. 1.3a) 
[ 2 ] 
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Fig . 1.3c. 



' S-v ■ •'.V **"' 




An example of two-dimensional filtering 
(High pass filtered version of Fig. 1.3a) 
[ 2 ] 
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nonrecursive one in number of computations and storage 
requirements. Hall demonstrates that the recursive filter 
algorithm constitutes the best method for large data, i.e., 
it is the fastest and cheapest. 

When processing two-dimensional data by recursive filters 
a fundamental problem exists due to inherent feedback, namely, 
the problem of numerical stability. Since, in general, two 
variable polynomials cannot be factored into a product of 
first and second order real coefficient polynomials in each 
of the variables, it is difficult to solve the stability 
problem for two-dimensional recursive filtering. Conse- 
quently, the majority of papers published in two-dimensional 
digital filtering deal with the design of nonrecursive 
filters which are inherently stable [6-10]. 

There are several papers discussing two-dimensional 
recursive digital filters. For example [11] and [12] 
formulate two-dimensional recursive filters by the z- 
transform and linear difference equations, and although 
they investigate problem areas related to stability and 
realization, the majority of problems remain to be solved. 

The most important applications of two-dimensional 
digital filtering of digital data are in picture processing 
and geophysical data analysis. Picture processing can be 
categorized into: 

a) Digital image restoration and enhancement, 

b) Computer pictorial pattern recognition. 
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The methods of image restoration and enhancement are 
applied to invert degradations, such as aberrations, 
atmospheric effect, scanning, motion, and to manipulate 
images to improve viewing phenomena experienced by the human 
eye. Image restoration and enhancement has been used with 
great success in biomedicine, i.e., in extraction of quan- 
titative information from x-ray films, chromosome counting, 
measuring the extent of arteriosclerosis from arteriograms, 
in forensic sciences application, i.e., fingerprint image 
enhancement for automatic classification, and in astronomy, 
i.e., removal of turbulence from astronomical photography 
[13] . 

Pictoral pattern recognition by digital computer has 
gained great importance in satellite surveillance and mili- 
tary reconnaissance. The application of two-dimensional 
digital filters permits the separation of different horizontal 
scans on magnetic and gravity maps in geophysics and can 
be used equally well for structural and topographic maps 
or for any other type of factor which is available in the 
format of a planar grid [14]. Although digital filtering 
in several dimensions is gaining importance in medicine, 
i.e., heart volume measurements, in fire control problems, 
and to analyze complex electronic circuitry, there exists 
no general theory concerning structures, analysis and 
design. 

These considerations can be summarized as follows. 

Since one-dimensional recursive digital analysis and design 
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techniques depend strongly on the fundamental theorem of 
algebra and on extensive utilization of continuous design 
theory/ it represents a special, i.e., a non-generalizable 
field in the area of N-dimensional recursive digital filtering. 
Multi-dimensional digital filtering is well developed in two- 
dimensions but is based, due to the absence of recursive theory 
analysis and design tools, predominantly on nonrecursive 
filtering schemes. Because of their inherent advantages, much 
of the current research is directed towards recursive filtering 
algorithms and techniques. 

B. AREA OF INVESTIGATION AND OBJECTIVES 

One area of investigation of this thesis is to see how the 
semi-infinite straight line approximation technique of one- 
dimensional filters (based upon the log modulus/log frequency 
or decibel vs. log-frequency plots) can be extended to two 
(or higher) dimensional recursive filters. The advantages of 
log-modulus approach is that a transfer function in factored 
form becomes a linear sum of terms when the logarithm of its 
magnitude is taken, and each of these terms can be approxi- 
mated by semi-infinite straight line segments in the one- 
dimensional case, and semi-infinite planes in the multi- 
dimensional case, when expressed in terms of log-frequency. 

The following are specific objectives: 

1. To postulate typical (canonic) factors for 
two-dimensional filters and to investigate their properties 
in the log-modulus/log frequency domain. 
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2. To develop an approximation technique for the 
design of two-dimensional recursive filters to meet given 
specifications in the frequency domain. 

The frequency response of two dimensional filters is 
inherently a volume in the space of magnitude, log and 
log This volume is the given or desired specification 

of the filter, which in the separable case can be approxi- 
mated by the combination of semi-infinite planes. 

In this thesis we investigate how we can approximate 
this volume by semi-infinite planes as an extension of 
the approximation used in one dimension, i.e., by semi- 
infinite lines. After determining the planar equations 
which perform the given specification of the filter, we 
investigate ways to get a transfer function; first, in the 
continuous frequency domain and second for two-dimensional 
recursive equations in the discrete variable domain. 

The results of this study will classify those charac- 
teristics that can be achieved by cascading simple canonic 
factors, and indicate limitations that are imposed by 
trying to achieve a particular frequency characteristic 
(such as the fan filter) in as simple form as possible. 

The approach taken here is to investigate the problem 
in the continuous domain first, and then study the trans- 
formation from continuous frequency (differential equation) 
domain to discrete time (sampled data) domain. 
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C. PREVIEW OF RESULTS 



We have proposed a geometric method to design two- 
dimensional recursive digital filters. This is a new 
method which is relatively simple and easy with respect to 
other design methods. We have investigated properties of 
some typical (canonic) factors and considered their stability 
properties. In Chapter II we go through the basic properties 
of two-dimensional filtering. In Chapter III we provide 
background in terms of planar geometry for the discussion 
in subsequent chapters. In Chapter IV we investigate the 
properties of typical factors and their stability in continu- 
ous variable domain. The new design technique, including 
stability considerations, starting with the continuous 
frequency domain and extending to two-dimensional recursive 
filters in the discrete domain, is presented in Chapter V. 

The results are summarized and some comments for further 
study are made in Chapter VI. It should be noted that the 
results presented apply to quarter plane filters, as defined 
in Chapter V. 
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II. FUNDAMENTALS OF TWO-DIMENSIONAL RECURSIVE FILTERING 



A. INTRODUCTION 

There are many signals that are inherently two-dimensional 
for which two-dimensional signal processing techniques are 
required. The most common examples are, of course, images 
in which the two variables are the spatial coordinates. 

Image processing [15], [16] plays an important role in many 

areas of scientific and technical research. Some examples 
are satellite imagery radiography, radar, and biomedical 
images such as acoustical holograms , medical x-rays , and 
electron micrographs [2] . Not all two-dimensional data 
come from an image. In prospecting for gas and oil, a 
common procedure is to monitor the seismic signals produced 
by detonating an explosive charge [17] . An array of geo- 
phones situated along a line radiating in a vertical or 
horizontal direction from the point of explosion is used 
to record these signals. In this case, the two variables 
are distance (from the point of explosion measured along 
the line) and time (the time at which the seismic signal 
reaches a given geophone or, equivalently, a given distance 
from the explosion) . 

Although two-dimensional signals may be processed by 
one-dimensional systems, it is generally preferable to con- 
sider using two-dimensional systems. - Many of the basic 
ideas of one-dimensional signal processing may be readily 
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extended to the two-dimensional case. There are some 
very important concepts of one-dimensional systems, however, 
that are not directly extentable to two dimensions. It 
is the goal of this chapter to discuss the basic ideas and 
techniques of two-dimensional signal processing and to 
illustrate them into the context of two-dimensional filter 
design. A more detailed discussion of this theory plus 
an overview of recent work in two-dimensional filtering is 
contained in the excellent survey paper by Mersereau and 
Dudgeon [18], and H. Chang and J.K. Aggarwal [19]. Other 
useful references are the first few chapters of [1] and 
chapter 7 of [2] . 

J 

B. DEFINITIONS 

1. Two-Dimensional Signals 

One-dimensional signals are functions of a single 
variable and two-dimensional signals are functions of two 
integer variables. Consider the two-dimensional sequence 
x(m,n) where m and n are integer variables, shown in figure 
2.1. As in the one-dimensional case, the notation x(m,n) 
is often short hand for a sampled version of a continuous 
two-dimensional signal x(s,t): i.e.. 



x(m,n) = x(mT^,nT 2 ) 



x(s ,t) 



s=mT^, t=nT 2 



( 2 . 1 ) 



One interpretation of x is a function which assigns a 
(generally complex) number to each integer ordered pair 
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Fig. 2.1. Pictorial representation for two- 
dimensional discrete array 
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(irt,n) . The function x is not defined when either or both 
of its arguments (m,n) is not an integer. Such a signal 
will be interchangeably referred to as either a two-dimensional 
array, a two-dimensional sequence or simply as a two-dimensional 
signal . 

Some useful two-dimensional sequences are defined 
below and shown in figures 2.2 to 2.4. These include: 

1. Digital Impulse or Unit Sample 



u Q (m,n) 



1 m = n = 0 

0 elsewhere 



( 2 . 2 ) 



2. Digital Step 



u_ 1 (m,n) 



1 m,n >_ 0 

0 elsewhere 



(2.3) 



3. Exponential 



x(m,n) = 



( m n 
I 3l a 2 



m,n >_ 0 
otherwise 



(2.4) 



4. Sinusoid (Complex) 



j (aa-,m+a) 2 n) 

x (m,n) = e < m,n < +°° (2.5) 



As seen above, the two-dimensional step is related to the 
two-dimensional impulse by 
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Fig. 2.2. Digital Impulse 
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Fig. 2.4. Exponential 



m 



n 



u -1 (m,n) 



l l u 0 ( m i' n i) 



(2.5a) 




2 . Two-Dimensional Systems 

A two-dimensional discrete system is defined as 
the unique transformation T that maps an input sequence 
(x(m,n)} and the initial condition sequence {s(m,n)} into 
an output sequence (y(m,n)} 



In signal processing applications, we assume zero initial 
conditions. With this assumption, the above relation may 
be shortened to 



y(m,n) = T[{x(m,n) } , {s (m,n) }] 



( 2 . 6 ) 



y ( m , n ) = T [ {x (m,n) } ] 



(2.7) 



The system characterized by T is 



(a) linear if 



T[{ax^(m,n) + bX 2 (m,n)}] = aT [ {x^ (m,n) } ] + bT [ {X 2 (m,n) } ] 



( 2 . 8 ) 



for arbitrary constants a and b. 



(b) shift invariant if 



y (irrm-^n-n^ = T [{x (m-m^,n-n^) }] 



(2.9) 



for all m^ and n-^. 
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A hardware and software system which is equivalent to the 
transformation T of (2.6) is called a digital filter. The 
class of complex-valued sequences which are square summable 
will play the role of input and output signals of the 
digital filter. Linearity and shift-invariance are indepen- 
dent properties of a system in that neither property implies 
the other. Thus the system 

T[x(m,n)] = h(m,n) x(m,n) 

is linear but not shift invariant and the system 

2 

T [x (m,n) ] = x (m,n) 

is shift invariant but not linear 

For two-dimensional discrete linear systems , the 
impulse response is defined as the response of the system 
at (m,n) to a two-dimensional digital impulse input at 
(m^,n^) with zero initial conditions, i.e., 

h(m,n;m^,n^) = T [ {u q (m-m^,n-n-]_) }] (2.10) 

where u Q (m,n) denotes the digital impulse. In particular, 
the impulse response of a two-dimensional discrete LSI 
(Linear Shift ^Invariant) system is 
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h (m,n;m 1 ,n 1 ) 



h (m-m^ , n-n^) 



or 



h (m,n) = T [ {u Q (m,n) } ] 



( 2 . 11 ) 



The response of a two-dimensional LSI* system with zero 
initial conditions is completely characterized by its 
impulse response. In fact using the identity 



00 

x(m,n) = £ x (m^,n^) u q ( m-m^ ,n— n^) 

m l' n l = -“ 

00 

= l x (m-m 1 ,n-n 1 ) u Q (m^^ ,n 1 ) 

m 1 ,n 1 =- 00 

we get 



y (m,n) = T[{ ’ l x (m 1 ,n 1 ) u q ( m-m 1 / n-n 1 ) }] 

m l ,n l = ”°° 

oo 

= £ x (m^ / n^ ) T [ { u Q (m-m^ , n-n^ ) } ] 

m 1 ,n 1 =-“ 

00 

= l x (m 1 ,n 1 ) h (m-m 1 ,n-n 1 ) 

m l ,n l =-C ° • 



*LSI = Linear Shift Invariance 
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oo 



y (m,n) 



l x(m-m 1 / n-n 1 ) h (m 1 ,n 1 ) 

m r n r" 



= x(m,n) * h(m,n) 



( 2 . 12 ) 



where denotes two-dimensional convolution. In other words, 

for LSI systems the basic convolution theorem is valid. 

3 . Causality, Separability, Stability 

In discussing various features of two-dimensional 
discrete systems, it is convenient to have spectral nota- 
tion [19] to represent certain regions of the spatial domain. 

1. Let S(a,S) denote a sector with the angular 
interval (a, B) • In polar form S(ct,8) may be described by 






S(a,S) = { (r , 0) | r >_ 0 , a < 0 < 8 } (2.13) 



where (a, 8) are measured in radians. In particular, we 
give the following notation for the sectors frequently 
encountered : 



S++ = S [ 0 , tt/2 ] ; S-+ = S[tt/2,tt] 

S — = S [-it , -tt/2 ] ; SH — = S[-tt/2,0] (2.14) 

S*+ = S [ 0 , tt) ; S*- = S [ tt , 0 ) 

where "[" and "]" imply including the boundary and ")" implies 
not including the boundary. 
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Notice that S++, S — h, S — , and S+- represent first, second, 
third, and fourth quadrants, respectively, and S*+ and S*- 
represent half-planes which are symmetric to each other with 
respect to the origin. (See fig. 2.5.) 

2. For a two-dimensional sequence (x(m,n)}, the 

set 



E x = { (m,n) | x (m, n) f 0} (2.15) 

is called the support of {x(m,n)}. 

3. Let (h(m,n)} be the impulse response of a two- 
dimensional discrete LSI system, and let E^ be the support 
of {h (m, n) } . 

Then, a two-dimensional LSI system is said to be 

causal ifE^ is the set S++, and semi-causal if E, is the 
h n 

set S* + . They are shown in Fig. 2.5. 

A two-dimensional discrete LSI system is said 
to be separable if its impulse response can be factored into 
a product of one-dimensional responses (Fig. 2.6); i.e., 

h(m,n) = h^ (m) (n) (2.16) 

If the equation (2.16) is not satisfied, the filter is said 
to be nonseparable . The advantages of separable filters is 
that two dimensional convolution (2.12) may be carried out 
as a sequence of one-dimensional convolutions. This can be 
seen by rewriting (2.16) as follows. 
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Fig. 2.5. Support Area of Impulse Response 

(a) Causal 

(b) Semi-causal S* + 

(c) Semi-causal S*_ 
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Fig . 2.6. 



Separable Digital Systems 
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oo 



y(m,n) = £ h 1 (m-^) h 2 (n^) x (m-m-^ , n-n-^) (2.17a) 

m^ , n^=-« 

OT CO 

= l h 1 (m 1 )[ l h 2 (n 1 ) x (m-ir^ , n— n 1 ) ] (2.17b) 

m^=-°° n^=-<» 

OO 

= l h 1 (m 1 ) a (m-m^ , n) (2.17c) 

m^=-°° 

where a(m-m-^,n) is a sequence of one-dimensional convolutions 
(obtained by evaluating the terms inside the bracket of 
(2.17b) for each fixed value of m-^) . Equation (2.17c) shows 
that y(m,n) may be obtained by a second sequence one-dimensional 
convolution. 

If both the input sequence x(m,n) and the filter 
impulse response h(m,n) are separable, then it is readily 
seen that the output sequence y(m,n) is also separable. 

In this case we obtain the result [using (2.12) and (2.16)] 

OO oo 

y(m,n) = [ \ h 1 (n^) b (m-m-^) J [ \ h 2 (n 1 ) c (n-n^ ] (2.18a) 

m^=-°° n^=-°° 

= a (m) 8 (n) (2 . 18b). 



where 



x (m, n) = b(m) c(n) 
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A two-dimensional LSI system is said to be (BIBO) 
stable if every bounded input sequence produces a bounded 
output sequence. This condition is satisfied if and only 
if the two-dimensional impulse response of the system is 
absolutely summable, i.e.. 



I |h (m,n) | < co (2.19) 

m, n=-°° 

As in the one-dimensional case, (2.19) can be shown to be a 
necessary and sufficient condition for stability [20] . One 
problem with (2.19) is that it can be quite difficult to 
evaluate for an arbitrary h(m,n). 

4 . Two-dimensional Difference Equations 

As in the one-dimensional case, two-dimensional LSI 
systems can often be described by a two-dimensional linear 
constant-coefficient difference equation relating the output 
y(m,n) of the filter to its input x(m,n). The most general 
form for such a difference equation is shown by 



l a (k 1 ,k 2 )x(m-k 1 ,n-k 2 ) 

(k l' k 2 )£R a 



l b(^ 1 ^ 2 )y (m-5, 1 ,n-il 2 ) 

(^l' ^ 2 ) eR b 



( 2 . 20 ) 



0 



where R^ and are finite sets of spatial grid points, 
called the input and output mask, respectively, and (a(k^,k 2 )} 
and {b(Z 1 ,Jl 2 )} are the set of constant coefficients that 
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characterize the particular filter. Generally, a difference 
equation has a family of solutions as is the case with a 
differential equation. The difference equation is said to 
be recursive if there exists a sequence of computations which 
yields the output sequence serially from the input sequence 
and that portion of the output sequence which is already 
computed, including the initial condition. An implicit 
ordering of grid points in the spatial domain is assumed for 
the serial computation of the output sequence. 

There are two cases of (2.20): 

1. R a , R^ in the S++ (first quadrant), b(0,0) ? 0 
and y(m,n) = 0 for outside of S++ (zero initial conditions) , 
then (2.20) represents a 'causal system. 

2. If R , R, in S*+, b(0,0) ^ 0 and y(m,n) = 0 for 

3. O 

outside of S*+ (zero initial conditions), then (2.20) 
represents a semi-causal system. 

In both cases we may assume b(0,0) = 1 without loss of 
generality and write (2.20) as 



y(m,n) = £ a (^ ,k 2 ) x (m-k^ ,n-k 2 ) 



(k 1 ,k 2 ) E R a 



l b(£ 1 ,A 2 )Y(m-£ 1 ,n-£ 2 ) 

(£ 1 ,£ 2 ) eR b -(0,0) 

( 2 . 21 ) 



It is clear that (2.21) is a recursive equation, and the 
implicit ordering, for example may be given as 
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{ ( 0 , 0 ) ; ( 1 , 0 ) , ( 0 , 1 ) , ( 1 , 1 ) ; ( 2 , 0 ) , ( 2 , 1 ) , ( 0 , 2 ) , ( 1 , 2 ) ; ( 2 , 2 ) . .} 



for a causal system, and 

{ ( 0 , 0 ) , ( 0 , 1 ) , ( 0 , 2 ) (- 1 , 1 ) , ( 0 , 1 ) , ( 1 , 1 ) . 

... (- 1 , 2 ) , ( 0 , 2 ) , ( 1 , 2 ) , . . . ; ...} 

for a semi-causal system. 

If {b(&^,$, 2 )} = {u o (i^,£ 2 )} the equation (2.21) 

becomes 

y(m,n) = £ a (^ , k 2 ) x (m-k.^ , n-k 2 ) (2.22) 

(k l' k 2 )eR a 

so that the output sequence is a weighted moving average of 
the input sequence. The impulse response of (2.22) is 
unique and is obviously {a(k^,k 2 )}. Such filters are called 
finite extent impulse filters (FIR) because the impulse 
response has only finitely many nonzero terms. FIR filters 
obviously are always BIBO stable because a(k^,k 2 ) is always 
absolutely summable. (It has only finitely many nonzero 
terms . ) 

If, on the other hand, the sequence {b(i^,^ 2 )} 
has several nonzero terms, then the impulse response will 
generally have infinitely many terms. We will say that such 
filters are infinite-extent impulse response filters (HR) . 
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Such filters may be BIBO unstable, and so it is necessary 
to devise a stability test for such filters. 

FIR filters have the advantages that they are 
easier to design, that controlling the phase response is 
easier (in particular, it is easy to produce filters which 
have a linear phase response) , and that since the input- 
output relationship is a finite-extent convolution. Fast 
Fourier Transform (FFT) techniques can be used to speed the X 

computation of the output sequence. On the other hand, 
recursive (HR) filters generally require less high speed 
storage and usually require less computation time than FIR 
filters. Comparison of FIR and IIR filters is contained 
in [5]. 

5 . Two-dimensional z-tranform 

A useful mathematical tool for representation of a 
sequence x(m,n) is its two-dimensional z-transform. Given 
a sequence x(m,n), the z-transform of this sequence is 
defined to be: 

< 

CO CO 

Z [ {x (m,n) } ] 4 X(z 1 ,z 2 ) = 11 x(m,n)z^Z 2 (2.23) 

m=-<» n=-°° 

where z^ and z 2 are independent complex variables. The most 
useful property of the z-transform is the fact that the z- 
transform of the convolution of two sequences is the product 
of their z-transforms , i.e.: 
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Z [ { x (m, n) }*{h(m,n) }] 



(2.24) 



= X(z 1 / z 2 )H(z 1 ,z 2 ) 

Therefore, the input-output relationship of an LSI filter 
may be expressed in the z-transform domain as: 

Y( Z;l ,z 2 ) = H(z 1 ,z 2 )X(z 1 ,z 2 ) (2.25) 

where H(z^,z 2 ), the z-transform of the impulse response 
is called the transfer function. 

A two-dimensional transfer function evaluated on 
the unit circle is called the frequency response of the 
system, i . e . , 



'i 

jw, ju) 9 jto. joj 9 jco. 

Y (e ,e ) = H (e ,e )X(e ,e ) (2.26) 

when a two-dimensional LSI system is describable by the 
difference equation of (2.20), and we take a z-transform 
of this, it follows that 

— k] — k„ 

[ l a(k 1 ,k 2 )z 1 "z 2 ]X(.z 1 ,z 2 ) 

(k 1 ,k 2 ) sR a 

-Z -Z 

- [ l bU v Z 2 )z 1 1 z 2 2 ]Y( Zl ,z 2 ) = 0 

(^l' eR b 
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or 



H(z r z 2 ) 



Y (z 1 , z 2 ) /X (z 1 , z 2 ) 



“ k i _k o 

I a(k 1 # k 2 )z 1 1 z 2 
( k l'k 2 ) £R a 



l b(£ 1 ,£ 2 )z 1 1 z 2 

d 1 ,l 2 ) ERj^ 



A(z 1 ,z 2 )/B(z 1 / z 2 ) 



(2.27) 



where A and B are the z-transf orms of {a(k^,k 2 )} and 

respectively. In particular when (2.27) repre- 
sents a causal or semi-causal system, i.e., (2.27) is the 

z-transform of (2.21), we can give explicit descriptions for 
R and R, in each case. For the causal case we define 



a 



b 



(Fig. 2.5a) 



R 



'a 



{ (k^ , k 2 ) | 0 <_ k-^ <_M^ , 0f. k 2 — 



(2.28) 




For the semi-causal case we define (Fig. 2.5b) 



R 



a 




(2.29) 



R, 



L b 
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Fig. 2.5b. Diagram of how output sequence is computed with 
semi-causal filter 



and in both cases b( 0 , 0 ) — 1 is assumed. The inverse 
z-transform is defined as 



x(m,n) = 5- / / X(z.,z ) z™ 1 z!J 1 dz.dz. ( 2 . 30 ) 

(2njT 12 12 12 

C 1 C 2 

where c-^ and C2 are suitable closed contours in the z-^ and 
Z2 planes. 

The two-dimensional Fourier transform of a sequence 
is defined as 

jw-, jtOp 

X (e , e ) = F [x (m,n) ] 



X (z 1 / z 2 ) 



( z L , z 2 ) sr‘ 



- j (aj. m+to~n) 

l x (m,n) e ( 2 . 32 ) 



m, n=-°° 



For any given sequence the set of (z^,Z2> for which the 
z-transform converges is called the region of convergence. 
Uniform convergence of the z-transform requires that the 
sequence be absolutely summable, i.e. 



I |x(m,n)||z 1 | m |z 2 | n < 00 ( 2 . 32 a) 

m, n=->» 
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3 3 

Similarly, the Fourier transform X(e ,e ) converges 
uniformly to a continuous function of and if 



l |x(m,n) | < co . (2.32) 

m, n=-°° 



A power series of the form of (2.23) is known as a two- 
dimensional Laurent series in the theory of complex functions 
in two variables . 

6 . Two-dimensional DFT 

Let (x(m,n)} be a finite support sequence of the 
size MxN. Then the two-dimensional discrete Fourier 
transform (DFT) of {x(m,n)} is defined as 



X(k x/ k 2 ) 



M-l N-l 

l l x (m,n) e 
m=0 n=0 



+ g* 2 > 



(2.33) 



for 0 <_ k 1 < M-l and 0 < k 2 < N-l. Similarly, the inverse 
discrete Fourier transform (IDFT) is given by 
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X(m,n) = 



MN 



M-l 

I 

V 0 



N-l 

l 

k 2 =0 



x(k 1 ,k 2 )e 



ki ^2 

+j2Ti(— m + — n) 



(2.34) 



Note that the DFT corresponds to sampling the Fourier trans 
form of M x N points, i.e.. 



X(k 1 ,k 2 ) 



3“i 



jo). 



= x (e x ,e *) k x k 2 

“r 2,T M-'“2* 2,t ir 



(2.35) 



j“l 3^ 2 

More generally, if X(e ,e ) is the Fourier transform 
of (x(m,n)} whose size is greater than MxN, then the 
IDFT of {X(k^,k 2 )} will result in aliasing. In other words 
the IDFT of (X(k^,k 2 )} will yield one period of (x(m,n)} 
where 



oo 

x(m,n) = l x (m+m^M,n+n^N) (2.36) 

nu ,n^=-<» 
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C. STRUCTURES OF TWO-DIMENSIONAL FILTERS 

1 . One-dimensional Digital Filter Realization 

One-dimensional digital filters, in general, can 

be shown in the z-transform domain as follows: 

/ 



H(z 1 ) 



c . a 
l b £ 



1=0 



M 



1 - l a m z 



-m 



m 1 



m=l 



y( Z; l ) 

X(Z X ) 



(2.37) 



or, in the time domain 



y (n) 



M L 

l a^ y(n-m) + l x(n-£) 

m=l 1=0 



(2.38) 



The realization implementing equation (2.38) directly, is 
known as the direct form 1 and shown in Figure 2.7. 

Equation (2.37) can be rewritten in a slightly 
different form by introducing a new variable, W(z ][ ), such 
that 



W(z L ) Y(z 1 ) 
H(z l } = X ( z x ) * W(z x ) 



(2.39) 



The corresponding set of difference equations consists of: 
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Fig. 2.7. Direct Form 



M 

w(n) = x (n) + l a m w(n-m) (2.40) 

m=l 

and 

L 

y(n) = l b z w(n-£) (2.41) 

1=0 



Equations (2.40) and (2.41) are realized in the 
direct form 2 as shown in Figure 2.8 which can be visualized 
as a cascade realization of two digital filters, realizing 
the denominator and numerator polynomials, respectively. 
Figure 2.8 can be redrawn as shown in Figure 2.9. The 
resulting realization is called direct form 3 or canonic 
form, since it represents a structure with the number of 
delays equal to the order. 

At this point, an additional realization will be 
introduced, the direct form 4, which combines the charac- 
teristics of the direct forms 1 and 2 having only two 
entries in each summer, which corresponds directly to 
hardware implementation and by realizing the numerator and 
denominator polynomials separately in cascaded form 
(Figure 2.10) . 

There are several other methods to realize a one- 
dimensional M order digital filter. For example, the 
cascade form of first and second order sections H(z^) , as 
shown in Figure 2.11, where 
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Fig. 2.8. Direct Form 2 
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L 



H (z 1 ) 



l 

i=0 




M 



1 



I 

m=0 



a 



m 



z 



-m 

1 




Fig. 2.9. Direct Form 3 
(Canonic) 
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H (z 
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Fig. 2.10. Direct Form 



k 

H(z l ) = a 0 11 k = ^y=- (2.42) 

i=l 

Each component H^(z^) can be realized in one of the above 
outlined forms. 

Partial fraction expansion methods applied to 
equation (2.37) lead to a parallel form realization of 
first and second order sections (Figure 2.12), i.e., 

k 

H(z x ) = c + l H i (z 1 ), k = [~i] (2.43) 

i=l 

Other forms of realization include hybrid structures, 
i.e., parallel-cascade forms, the transpose configuration, 
which can be obtained for all previously outlined structures 
by reversing the direction of signal flow and by interchanging 
all branch and summing modes, wave-digital and continued 
fraction expansion filters. 

The direct forms discussed above and the series/parallel 
arrangement of lower order sections are by far the most 
often used ones in computer simulation and in digital 
hardware . 

2 . Two-dimensional Digital Filter Realization 

Two-dimensional quarter plane digital filters, in 
general form, can be shown as 
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Fig. 2.12. Parallel Form 



H(z 1 ,z 2 ) 




L 



2 




Z l l 2 





1 



M 



1 



M 



2 



l l 

m^=0 m 2 =0 



m l m 2 





(2.44) 



The numerator of the two-dimensional transfer function (2.44) 
can be written as a one-variable polynomial in z 2 , weighted 
by coefficients which are functions in z^, i.e., 



L 2 L 1 -% 2 
N(z x ,z 2 ) = l ( l z z i )z 2 (2.45) 

z 2 =° IfO 1 2 



2 

l 

i 2 =0 



-l. 



(Ni 1 ^2 (Z;l) } Z 2 



(2.46) 



The same reasoning applied to the denominator polynomials 
leads to 



D ( z ]_ ' z 2 ) 



M M 

1 1 -m. -m_ 

1 - l ' 2 a m^m 2 Z 1 )z 2 

m 2 =0 



(2.47) 



M, 



1 - 



l ^ D M.m_^ Z l^ Z 2 



-m 



m 2 =0 



M l m 2 



(2.48) 
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To obtain the direct form 3 representation, a new function, 
W(z 1 z 2 ) is introduced: 



H < z l z 2> 



^ ( Z 1 ' z 2 ) Y(z^,z 2 ) 

X(Zi,z 2 ) W(z-^ , z 2 ) 



(2.49) 



where each factor can be written using equation (2.46) as 



W(z r z 2 ) 



i 2 =0 



l <Ny 2 <V>* 2 2 X(z l' z 2 : 



(2.50) 



and equation (2.48) as 



M, 



-m. 



Y ( z i ' z 2 ) ^ (Z]_ ' z 2 ) £ (Dj„j ^ (z ^)) z 2 Y ( z ^ z 2 ) 

— n 1 ^ 



m 2 =0 



(2.51) 



The realization of equations (2.50) and (2.51) is shown in 
Figure 2.13 using compact notation, and Figure 2.14 where 
N l £ (z x ) and 

° M l m 2 ( z^) are implemented for each (L^S^ 
and (m^m 2 ) by one-dimensional direct form 4 realizations. 

The general form of two-dimensional causal digital 
filter (2.44) can be rewritten by introducing a new variable, 
W(z x z 2 ) or in the previous function, then, 



W(z 1 z 2 ) Y(z 1 ,z 2 ) 

H(z 1 /Z 2 } = X(z x ,z 2 ) * W(z 1 ,z 2 ) 



(2.52) 
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Fiy. 2.13. Direct Form 3, 2D. (Compact Notation) 
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Direct Form 



Each factor can be written as: 



Y(z 1' z 2 ) 



L 1 L 2 

I l b 

Al=° i 2 =0 



^1 ^2 

£ 1 £ 2 Z 1 Z 2 W(z l' z 2 } 



and 



(2.53) 



M, 



W(Z 1 /Z 2 ) =,r X(Z 1 ,Z 2 ) + l 

m i = 

(m 1 +m 2 ) t*0) 



M 
l a. 



mi=0 m 2 =0 



-m, -iru 

_ _ Z, Z W ( z. , Z_ ) 

m l m 2 12 12 



(2.54) 



The equations (2.53) and (2.54) can be written as a differ- 
ence equation as follows. 



y (m,n) 



L 



2 





b^ % W(m-5, 1 ,n-il 2 ) 



(2.55) 



and 



w(m,n) = x(m,n) + \ ^ a m m w (m-m 2 / n-m 2 ) (2.56) 

m^=0 m 2 = ^ ^ 2 

m l ,m 2 ^ o 
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Without loss of generality, is chosen equal to ]XL for 
all i. The structure realizing equation (2.55) and (2.56) 
is shown in Figure 2.15. It is noted that the number of 
delays in Figure 2.15 equals the order of H(z 1 ,z 2 ). The 
structure is, therefore, by definition, canonic. 

The direct form 4 (Fig. 2.10) realization of H(z^z 2 ) 

is shown in Figures 2.16 and 2.17, using compact notation and 

one-dimensional direct form 4 (Fig. 2.10) implementation of 

each N . and D M , respectively. 
l 1^2 M l m 2 

D. STABILITY 

Two major problems in the design of a recursive filter 
are approximation and stability. Since the output values 
are used through feedback by the recursive algorithm, it is 
possible for the output values to become arbitrarily large 
independent of the input. A filter of this kind is said to 
be unstable, which is an undesirable condition. Thus we 
need to know what constraints to put on the recursive filter 
coefficients, so 'that the filtering operation will be stable. 
This problem is similar to the stability problem for the one- 
dimensional case, except that the added dimension inherently 
increases the complexity of the analysis. The major diffi- 
culty is that the concept of zeros and poles does not hold 
in multi-dimensional systems so that simple algebraic systems 
extrapolation of one-dimensional results is not obvious. 

! Stability Theorems : 

Theorem 1 (Shank's theorem): A causal recursive filter 

with the z-transform H(z^,z 2 ) = A (z^ , z 2 ) /B ( z^ , z 2 ) , where A 
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Fig. 2.15. Direct Form 3, 2D (Canonic) 
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Fig. 2.16. Direct Form 4, 2D (Compact notation) 
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Fig. 2.17. Direct Form 







0 



and B are polynomials in z 1 and z 2 , is stable if and only 
if there are no values of and z 2 such that B(z 1 ,z 2 ) = 
for all z^ > 1 and z 2 _> 1. (In general we assume 

B(z 1' z 2 ) = b 00 + b 10 Z l + b 01 Z 2 + b ll z l lz 2 1 

that is an expression in increasing powers of z^ and z 2 .) 

In other words, the theorem says that if there are any 
values (real or complex) of z ^ and z 2 for which B(z^,z 2 ) 
is zero for z^ and z 2 simultaneously greater than or equal 
to unity in magnitude, then the filter will be unstable. 

It is much more difficult to determine the stability of 
two-dimensional filters than of one-dimensional filters. 

In the one-dimensional case, it is only necessary to locate 
finite set of roots in the z-plane. In order to mechanize 
the theorem, we must, in general, find the values of z^ and 
z 2 for which B(z^,z 2 ) = 0. There is no technique for locating 
the zeros of a general two-dimensional polynomial in (z^,z 2 ). 
Shanks, Tretial and Justice, in their paper [11] proposed a 
technique to apply theorem 1 to determine the stability of 
a two-dimensional recursive filter, H(z^,z 2 ), by which the 
unit disk in the z^ plane must be mapped into the z 2 plane 
by solving the implicit two-variable denominator polynomial 
of H(z^,z 2 ). A necessary and sufficient condition for 
stability is determined if the map of the unit disk 
di = (z x , z 2 _> 1) does not overlap the unit disk 
d 2 = (z 2 ,z 2 - in the z 2 P lane * This method requires an 
infinite number of mappings and, therefore, can not be 
applied exactly. 
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Huang [20] simplified Shank's stability theorem 
considerably by stating the stability theorem as follows: 
Theorem 2: A causal filter with a z-transform 



H ( z ! ' z 2 ) 



A(z 1 ,z 2 ) 
B ( z x , z 2 ) 



where A and B are two-variable polynomials in z^,z 2 , is 
stable if and only if: 



1) the map of the unit circle 6d. 



(z 1 , | z ± |= 1 ) , 



according to B(z^,z 2 ) = 0, lies outside of the unit disk 

^2 ~ ( z 2 ' I 2 2 I — ^ ; ’‘ n the z 2~P-*- ane ' anc ^ 

2) no point in the unit disk d^ 5 (z^,|z^|_>l) 
maps into the point z 2 = 0 by the relation B(z^,z 2 ) = 0. 

In order to apply this theorem as a stability test, we 
have to map the unit disk 6d^ = (z^,|z^|=l) into the z 2 
plane according to B(z^,z 2 ) = 0 or z 2 = f(z^) | | z l=1 and 
to see whether the contour in the z 2 plane lies inside the 
unit disk d 2 = (z 2 ,|z 2 |<l). Also, it is necessary to solve 



B(z^,0) = 0 to see whether there are any roots with magnitude 
greater than 1. 

Although theorem 2 is much simpler than using Shank's 
original theorem, the procedure is still infinite in the 
sense that, in principle, we have to map the unit circle 
<5d^ into the z 2 ~plane. 

Ansell [21] , whose main contribution is to couple the use 
of a Hermite test for checking stability with a series of 
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Sturm tests checking positivity, has reduced theorem 2. 
Although An sell's results enable us to test stabilty in a 
finite number of steps, it still is, unfortunately, very 
tedious . 

Make the change of variables: 



and. 



?2 



1 - 
T~+ 



z 



z 



2 

2 



and let 



H(z 1 ,z 2 ) 



E (P 1 /P 2 ) 
F (P 1 ,P 2 ) 



where E and F are polynomials in p^ and p 2 - We can restate 
Theorem 2 as follows. 

Theorem 3 (Ansell's theorem): The causal recursive 

filter H(z^,z 2 ) is stable if and only if: 

1) for all real finite to, the complex polynomial 

in p 2 , F(joj,p 2 ) has no zeros in ReP 2 > 0, and 

2) the real polynomial in p^, F(p^,l) has no 

zeros in Rep^ 0. 

This theorem is essentially the same as Theorem 2, but 
an advantage of Theorem 3 is that condition 1) can be 
tested using standard techniques of circuit theory. 
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Anderson and Jury [22] modified Huang's method by out- 
lining a procedure which replaces the bilinear transforma- 
tions by the construction of a Schur-Cohn matrix and checking 
for positivity of a set of self-inversive polynomials. It 
also replaces the Hermite test by a Schur-Cohn matrix test 
and requires a series of Sturm tests or equivalently, a 
system of tests establishing the roots distribution of a 
polynomial. These modifications represent a substantial 
reduction in computations as compared to Huang's method. We 
will discuss stability considerations which are related to 
our design technique in Chapter IV and Chapter V in detail. 

In this chapter, basic properties of two-dimensional 
digital filters have been reviewed and summarized including 
semi-causality to establish a fundamental theory of two- 
dimensional digital filters. 
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III. MATHEMATICAL BACKGROUND FROM ANALYTIC PLANAR GEOMETRY 



In this study, we propose using semi-infinite planes to 
approximate a given frequency response of a two-dimensional 
analog filters by extending one-dimensional semi-infinite 
straight line approximation to two-dimensional cases. For 
example, as we will see, the frequency response (dB vs. log- 
frequency) of the separable canonic factor 

H ( j U3 1 , jco 2 ) = (l+ju 2 x 2 ) 

can be approximated by a collection of planes. This is 
discussed in Chapter IV in detail. 

To prepare for the discussions in the next chapters, 
it is worthwhile to review and develop useful equations 
from analytic planar geometry which are going to be used 
throughout this study with some examples of how to use them. 

The properties of the analytic planar geometry can be 
summarized as follows: 

1. The general equations for a plane are given by 

Ax+By+Cz+D=0 (3.1) 

where A, B and C are the direction numbers of the normal 
to the plane. The equation of a plane passing through the 
point (x^,y^/Z^) and perpendicular to the line which has 
direction numbers of A, B and C is given by: 
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0 



(3.2) 



A(x-x 1 ) + B(y-y 1 ) + C(z-z 1 ) = 

2. The direction cosines of the normal to the plane 
are given by 



cos a 



cos 6 



cos y 



A 



(A 



2 + b 2 + c 2 ) 1 / 2 



(a 2 + b 2 B + c 2 ) 1/2 



(A 2 +B 2 + C 2 ) 1/2 



(3.3a) 



(3.3b) 



where a, 3, and y are the angles that the normal makes with 
the x, y, z axes respectively. 

The relationship between direction cosines is given by 



2 2 2 
Cos a + Cos 3 + Cos y = 1 



(3.3d) 



3. For a line joining points (x^,y^,z^) and ( X 2 >Y 2 ' Z 2 ) 
the equation is given by 

Cos a _ Cos 3 _ Cos y (3 4 ) 

x 2 - X 1 y2‘ y l z 2 ” Z 1 

4. The equation of a straight line with the direction 
numbers a, b, c and passing through the point (x^,y^,z^) is 
given by 



y-y x 

b 



2-2 



(3.5) 



X - 



a 



1 
c 



5 . The angle 0 , between two intersecting lines with 
direction numbers a, b, c and a', b', c' respectively is 
given by: 



Cos 0 = Cos a Cos a ' + Cos 6 Cos 0 ' + Cos y Cos y ' 

(3.6a) 



or 



Cos 0 



aa 1 + bb ' + cc 1 

Va 2 + b 2 + c 2 Va ' 2 +b' 2 + c'2 



(3.6b) 



6. For parallel lines, the relationship between 
direction numbers is given by 




(3.7) 



7. For perpendicular lines: 

aa' + bb' + cc' = 0 (3.8) 

8. The direction numbers of the line of the intersection 
between two given intersecting planes, with direction numbers 
A, B, C and A', B', C' respectively, is given by 
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a = 



B 



B' C' 



(3.9a) 



b = 



(3.9b) 



c = 



C' A’ 



B 



A’ B' 



(3.9c) 



The angle 9 between these two planes is given by 



Cos 9 



AA' + BB' + CC' 

'\Ja 2 + b 2 + c 2 V A ' 2 + B ' 2 + c ' 2 



(3.10) 



9. The equation for a plane in terms of its slope can 
be calculated as follows : 

As seen from Fig. 3.1, the slope, m, of the normal 
line is given by 



o h C 

m = tan 9 = — = 



(3.11) 



Va 2 ' ~ 2 



I A + B 



m = slope of the plane 



M = tan (90° - 9) = 



_ 2 = _ I 

h m 



(3.12) 
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z 




Fig. 3.1. A General Plane Configuration 



73 



The equation of the plane is given by one of the 



following : 




A(x-x^) +By + Cz = 0 


(3.13a) 


Ax + B(y-y^) + Cz = 0 


(3.13b) 


Ax+By + CCz-z^) = 0 

and 


(3.13c) 


Ax^ = By^ = Cz^ = -D 


(3.13d) 



The intersection of the plane with the z-x plane 



is given by 




A (x - x^) + Cz = 0 

and the slope, m^, of this line is 


(3.14) 


dz A 

m l dx C 


(3.14a) 



Similarly, the intersection of the plane with the 
z-y plane is given by 

B (y - y 1 ) + Cz = 0 (3.15) 
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and the slope of this line, m 2 , 



is 



IT^ 



dz _ B 
dy C 



(3.15a) 



Substituting 
yields the following 



(3.15a) and (3.14a) into (3.13) 
alternate equations for the plane. 



-m(x-x^) - ra 2 y + z = 0 



-m^x - m 2 (y-y 2 ) + z = 0 



-m^x - m 2 y + (z - z^) = 0 



(3.16a) 

(3.16b) 



(3.16c) 



The direction cosines of the normal in terms of 



m^ m 2 are 



-m. 



Cos a = 



Vm^ 2 + m 2 2 + 1 

-nu 



Cos 3 = 



Vm^ 



2 + m 2 + 1 



(3.17a) 



(3.17b) 



Cos Y 



Vm^ 2 + m 2 2 +1 



(3.17c) 



Substituting (3.14a) and (3.15a) into (3.11) and 
(3.12) yields 



m = 



(3.18a) 






2 2 
+ m 2 
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and 



M 



= - 



2 , 2 
+ m 2 



C3.18b) 



The angle, 0, between two intersecting planes from 
(3.10), (3.14a) and (3.15a), is given by 



Cos 0 - 



m 1 m 1 ' + m 2 m 2 ' + 1 



(3.19) 



V^ 2 + m 2 2 + 1 V m { 2 + m 2 2 + 1 



The direction numbers of the line of intersection, 
from (3.9), (3.14a) and (3.15a) are given by 



a = (m 2 - m 2 ) CC' 



(3.20a) 



b = (m£ - m^) CC' 



(3.20b) 



c = (m^m 2 - m 2 m^) CC' 



(3.20c) 



The direction cosines of this intersection line 
are given by 



Cos a = 



m 2 • m 2 



(3.21a) 



V 2 2 

(m£ - m 2 ) + (m£ - n^) + (m^m 2 - m^m^) 



Cos 3 = 



m i - m l 



V 2 2 2 

(m 2 -m 2 ) + (m£ -n^) + (m|m 2 -m^n^) 



(3.21b) 



76 



Cos y 



(3.21c) 



m l m 2 - m 2 m i 



V 



2 2 2 
(m^ -n^) + + ( m { m 2 



10. When two planes are added in one dimension to form 
a new plane, the following results. Given two planes 



-m. x - rruy + z - z, = 0 



-m^x - m^y + z' - z^ = 0 



Their sum in the z-direction with z" = z' + z is 



given by 



-(m^+m')x - (i^+n^'Jy + z" - (z-^ + zj^) - 0 



The slopes of the intercepts with the zx and zy planes 



m. 



= - (m 1 + mp 



m. 






(m 2 + m 2 ' ) 



The slope of the normal is given by 



m" = ( (m x + m|) 2 + (m 2 +m£) fc ] 



... 2 .- 1/2 



(3.22) 
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To aid in comprehension, two examples follow: 

Example (1) : Given m^ (slope in z,x plane) , m^ 

(slope in x,y plane) , and x^ (the point at which x-axis goes 
through the plane). Find the equation of the plane. 

The general plane equation (3.1) is given by 

Ax + By + Cz + D = 0 

When z = 0 



X = 


-i*-| 


m 3 “ 

and 


f - -! (given) 


x i = 

When y = 0 


- ^ (given) 


z = 


-§*-§ 


m l = 

and 


t§ - -i <^ en > 




D 



Z 1 C 
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If we substitute these into the general equation, the result is 



Ax - Am^y - 



m. 



-z - Ax, 



0 



or 



x - m^y 




x 



1 



0 



or 



z = m-^x - m^m^y - x^m^ (3.23) 

As shown here, we can derive the plane equation easily 
from the given specification. After finding this equation, 
we can find the canonic factor by simply letting 

x = log 

y = log ^2 

z = 10 log | T| 2 



Alternately, we can go from the canonic factor to equations. 
Thus equation (3.23) can be written as 



10 log | T ( jio,, j tu 2 ) 



where 



m^ log “ m^m^ log w ^ ~ 10 log K 

(3.24) 



10 log K = 
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This result assumes that 



|T(jw 1 / ju 2 ) 




(3.24a) 



Taking the logarithm of (3.24a) yields 
2 

10 log j T | = 10 p log id^ - 10 q log co 2 “ 10 log K 

(3.25) 

By equating (3.25) and (3.24) , we get 



(3.26a) 
(3.26b) 

It follows from (3.26a) and (3.26b) that 



10 p = m^ 

10 q = m^m^ 



= q/p • 

Thus the transfer function of (3.24a) has a logarithmic 
characteristic as sketched in Fig. 3.2. 

The previous result can be extended to the following 
transfer function 



T ( j / j co 2 ) 



(1 + jco, ) 



P/2 



-1/2.. . • . q/2 

K (1 + Ju> 2 ) 



(3.27) 



By substituting (3.26a) and (3.26c) into (3.27) , then 
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a 



Fig . 





3.2. The logarithmic characteristic of 

P 



I T ( ju 1 , j to 2 ) 



U). 



Kco. 



20 log | T | = m 1 log log co 2 - 10 log K 
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T(jw 1 ,j(o 2 ) = 



(1 + jw^) 



m /20 



K 1/2 (l + j oj 0 ) 



m 2 /2 0 



When oj 1 << 1 and w 2 << 1 , region I 



T(jw 1 / joj 2 ) = 



= K 



- 1/2 



and 



20 log | T | = - 10 log K 



When u, << 1 and oj 2 >> 1, region II 



T ( j / j to 2 ) 



1/2 . . ™ 2 /20 
K (3« 2 ) 



and 



20 log | T | = - 10 log K - m 2 log gj 2 



When u> 2 << 1 and >> 1, region III 

m, /20 

(jw x ) 

T ( j w ]L / j u> 2 ) = ~T/2 

K 

and 

20 log | T | = - 10 log K + log o> 1 



(3.28) 
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When >> 1 and w 2 >> 1, region IV 



T(j(D 1 ,ju) 2 ) 



Du 



m^^/20 



1/2 m 2 /20 

K (Ju 2 ) 



and 



20 log | T | = - 10 logK + m^log - m 2 log oj 2 . 



The resulting planar approximation for the logarithmic 
transfer function is shown in Figure 3.3 
Example 2 ; 

In this example, we derive planar equations from the 
given transfer function. 

Given a separable two-dimensional transfer function: 



T( joi,, ju 2 ) = ^ ^ (3 * 29) 

(1 + jcD^T^Ml + ja) 2 T 2 ) 4 

The magnitude square of the transfer function (3.29) is 
given by 



j T ( j a) ^ , j a) 2 ) | 



K 



,, A 2 2,p.. 2 2 <3 

(1 + co^ x ) ( 1 00 2 X 2 ) 



(3.30) 



Take the logarithm of both sides of (3.30), 



10 log I T I 2 = -lOp log (1 +(jo 1 2 t 1 2 ) - 10 q log (1 +(jj 2 2 t 2 2 ) + 20 log K 

(3.31) 



when u 1 t 1 << 1 and oo 2 x 2 << 1, the equation (3.31) becomes 
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*rtol<|lT| 





W 



t 




20 log |T | = -10 iLog K + m-^log - m 2 log w ^ 
84 



10 log I T I 2 = 20 log K = C 

where C is some constant. 

When oo^t^ >:> 1 and W 2 T 2 >> ^ ' 

10 log | T [ 2 = -20p log co^ - 20q log co 2 t 2 

+ 20 log K 



or 

2 

10 log | T | = -20p log co^ - 20q log co 2 “ 20p log 

-20q log x 2 + 20 log K 

Let 

z = 10 log | T| 2 



x - log 
y = log co 2 

and we get the planar equation, 



Z = -20px - 20qy - 20p log x^ - 20q log x 2 + 20 log K 



or 
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(~20p) x + (-20q)y - z 



or 



= 20 log 




q 



-) 



+ m 2 y - z = k 1 

where 

m 1 = slope of the line (x,z) plane 
= -20 p db/decade. 

m 2 = slope of the line (y,z) plane 
= -20q db/decade. 

In a similar manner the regions << 1; 

^ 2 t 2 >> ^ an< ^ ^“l 1 ! >:> U 2 X 2 <<: ^ can investigated. 

The planar approximation for the logarithmic transfer function 
is sketched in Fig. 4.1 where q = p = 1. 

In this chapter we have developed some useful planar 
geometry formulas and shown how they can be used to determine 
a planar approximation for the two-dimensional logarithmic 
transfer function for separable cases. Some non-separable 
factors are discussed in the next chapter in which experimental 
verification of this approximation technique is also presented. 
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IV. CANONIC FACTORS IN CONTINUOUS FREQUENCY DOMAIN 



In Chapter III, we discussed some equations of planar geo- 
metry and showed how to use them for approximating separable two- 
dimensional transfer functions. In order to improve the scope 
of the proposed design technique, we now consider the behavior 
of other separable and non-separable factors insofar as 
their planar geometry approximation is concerned. This will 
provide insight into the frequency response characteristics 
of two-dimensional systems and indicate what can be expected 
by cascading these canonic factors, knowing their individual 
behavior. The logarithmic characteristic of cascaded factor 
is simply the sum of their individual logarithmic character- 
istics. This approach parallels the semi-infinite line 
approximation (dB vs. log to ) technique which has been used 
so successfully in one-dimensional design. 

A. SEPARABLE FACTORS 

The following several separable canonic factors are now 
considered. 



( 1 ) 



(1 + j(0 1 T 1 ) P (l + jt 2 t0 2 ) q 



(4 .a) 



( 2 ) 




(3) 



[1 + (j(o 1 T 1 ) P ] a 



( 4 . c) 



where p, q and a are integers. 
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(1) First, as an example of the technique, consider 
the following analog separable transfer function, the first 
separable canonic factor with p = q = -1; 



H ( j CO , j 00 — ) = -TT ; T 7 - 7 -. : ; 7 r- 

1 2 (1 + (1 + joj 2 t 2 ) 



(4.1) 



The logarithm of magnitude-square of this function 
is given by 



10 log |h(co 1 ,oj 2 ) | 2 = “10 log (1 + - 10 log (1 + w 2 x 2 2 ) 

(4.2) 

Let us define the following for simplicity of notation. 



z = 10 log |h(oj^,co 2 ) 1^ = 20 log |H(co^,u) 2 ) | (dB) 

(4.3a) 

x = 20 log oj^ (4.3b) 



y = 20 log oo 2 (4.3c) 

Consider the regions in the (cjj^,oj 2 ) plane (as seen 
in Figure 4.2) as follows: 

Region I, when < 1 and to 2 x 2 < ^ ' 

then z=0 db (4.4a) 



Region II, when w^ T l > ^ an< ^ W 2 X 2 < ^ 

2 = -20 log (oo^t^) db (4.4b) 
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(4.4c) 



Region III, when < 1 and ^ 2 T 2 > 1 

z = -20 logtc^x-,) db 

and 

Region IV, when > 1 and ^ 2 X 2 > 1 

z = -20 log (o) 1 x 1 ) - 20 log (a) 2 x 2 ) (4.4d) 

The above equations (4.4a), (4.4b), (4.4c) and (4.4d) 

can be expressed as 



Region 


I , 


Z = 


0 






(4.5a) 


Region 


II, 


Z = 


-x - 20 


log 


T 1 


(4.5b) 


Region 


HI, 


Z = 


-y - 20 


log 


T 2 


(4.5c) 


Region 


IV, 


z = 


-x - y - 


20 


iog t 2 


- 20 log x -1 



(4 . 5d) 



These equations describe four intersecting planes in 
x,y,z space as shown in Figure 4.1. 

Plane I is the horizontal plane at zero db. Plane II 
and Plane III have slopes of -20 db/decade passing through 
the corner lines as indicated in Figure 4.1. These lines 
are analogous to the break point frequencies of one-dimensional 
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Fig. 4.1. Planar Approximation for 

H(aJ l ,aJ 2 } = ( 1+ jaj 1 T 1 )* (l+ja) 2 T 2 ) 

20 log | H (oo^ , co^) | = -10 log (l+w^x^) 

-10 log (I+oj^^t 2 ^ 
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°g w, 



Fig. 4.2. The regions of 

H(jw-j/jw 2 ) “ (1+jw^T^) ( 1+ j u 2 t 2 ) 

2 2 

20 log | H (co 1 ,o> 2 ) I = "10 log (1+e^ x ± ) 

-10 log {l+u^x^) 
in (log w./log u^) plane 
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log-modulus plots. Several properties of this plane are 
summarized. Pertinent equations from analytic geometry are 
given in Chapter III. 

The actual frequency response of this separable 
canonic factor with = T 2 = 1 is shown in Figure 4.3a 
(dB vs. log log ) which fits the planar approximation. 

In Fig. 4.3b, the contour plot of the frequency response is 
given, which verifies the regions in log cu^, log plane 
shown in Fig. 4.2. The logarithmic plots versus linear 
frequency are given in Fig. 4.3c and Fig. 4.3d. It is 
apparent that the approximation regions are well defined 
when logarithmic frequency is used. 

The frequency response of this separable canonic 
factor with = 1, T 2 = 2 is given in Figures 4.4. These 
results again confirm the logarithmic approximation. 

(2) Next consider the canonic form 

H(joi 1 ,ju 2 ) = (1 + jco^) (1 + ju> 2 T 2 ) (1 + ja> 2 T 3 ) 

where we have two separable factors of The frequency 

response of this factor with = 1, r 2 = ^ an< ^ T 3 = 4 is 
shown in Figures 4.5. As seen from this figure, we have 
-40 dB/decade slope in the direction of as expected 

from the one-dimensional case. The regions are shown in 
the contour plot with dotted line (Fig. 4.5c) . 
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In a similar manner we can analyze these regions as 
follows : The logarithm of magnitude square of this function 

is given by 

10 log |H(a) 1 ,u 2 ) | 2 = - 10 log(l + o) 1 2 t 1 2 ) -10 log (1 + u 2 2 t 2 2 ) 

- 10 log (1 + w 2 2 t 3 2 ) . 

Consider the following regions in (u)^,cd 2 ) plane (as 
seen in Fig. 4.5a) as follows: 

Region I, when << 1 and t0 2 T 2 < W 2 T 3 <<: ^ 

then z = 0 dB. 

Region II, when >> 1 and oj 2 t 2 < W 2 T 3 << ^ 

z = - 20 log dB 

Region II, when << 1 and ^2 T 2 < a) 3 T 3 <<: ^ 

z = -20 log (w 2 x 2 ) - 20 log(a) 2 T 3 ) dB 

Region IV, when >> 1, “2 T 3 > a) 2 T 2 >:> ^ 

z = -20 log w^x^ ~ 20 log ^2 T 2 ~ 2 ^* W 2 T 3 * 
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V, when <<: 1/ ^2^3 >> ^ 2^3 >>aj 2 T 2 >> 

Z = -20 logCo^T-j) - 20 logCo^tj) 

Region VI, when oo^t^ >:> W 2 T 3 > x 2 t0 2 >:> ^ 

Z = ~ 20 log ~ 20 log <*> 2^2 ~ 1°9 C0 3 T 3 

where 

Z = 20 log | H(w^/ t 2 ) I • 

(3) Consider the separable canonic factor now, 

H ( / ^2 ) = [ 1 + ( j T ]_ ) ^ ] 0 

We can show that the following regions occur by 
using the notation of (4.3a) and (4.3b) for simplicity: 
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z 






The frequency response of 

. _ 1 

H(jw 1 ,Du 2 ) (l+joi 1 ) (l+jo) 2 ) 

2Q log|H((o 1 ,a) 2 ) | = -10 log(l+u) 1 2 ) -10 log (l+u 2 ) 

(dB vs. log w^/log w 2 ) 



Fig. 4.3a. 



z = 
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I Oo. 



• s 



o 







QJ 



I. 



10 . 



I0O. 



X -SCALE= 1 . OOE+OO UNITS INCH. 
T-SCALE=1 . OOE+OO UNITS INCH. 



Fig. 4.3b. The contour plot of H((jj^,oj 2 ) = ' ( ~ i+ " jV 
z = 20 log j H | = -10 log (l+oo^ 2 ) - 10 



-► 

I og w 



1 

' (l + jcj 2 ) 

log (1 +oj 2 2 ) 
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z 







Fig. 4.3c. 



The frequency response of 
H(ju 1 ,ju 2 ) = {1+ja) .J (i+jo^) 

z = 20 log | H | = -10 logd+ojj^ 2 ) -10 log(l+w 2 2 ) 
(dB vs. (linear frequency) 
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Fig. 4.3d. 



The contour plot of H(w^,w 2 ) = 

z = 20 log | H | = -10 logd+u^ 2 ) 
(linear frequency) 



1 

(l+jw 1 ) ( 1+ j oj 2 ) 

- 10 log(l+oa 2 2 ) 
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Fig. 4.4a. The frequency response of 

H(ju 1 , jai 2 ) - (i+joj^) (1+j 2 u 2 ) 

2 2 
z = 20 log I H | = -10 logd+to. ) - 10 log(l+4w 2 ) 

(dB vs. log co^,log u^) 
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DO l " m bo 




X-SCRLE= 1 . 00E+00 UNITS INCH. 
Y-SCflLE= 1 . 00E+00 UNITS INCH. 



Fig. 4.4b. The contour plot of the frequency response of 
H(aj i t0 2 ) = (ir^qr d+j2cj 2 ) 

z = 20 log | H | = -10 log ( 1 +oj 1 2 ) - 10 log(l+4w 2 2 ) 
in (log o)^,log a> 2 ) plane 
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z 





Fig. 4.4c. The frequency response of 

H(ja) 1 ,jw 2 ) = ( i + j Ui ) 1 (i+j2u) 2 ) 

z = 10 log | H | = -10 logd+u^ 2 ) - 10 log(l+4oj 2 2 ) 
(dB vs. w^,a) 2 ) (linear frequency ccale) 
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>48.0 




Fig. 4.4d. The contour plot of the frequency response o^ 

. . 1 

H ( joo. , jw.,) = 5 T~ 

1 2 (l+co 1 2 ) (1+4 o) 2 2 ) 

z = 20 log | H | = -10 log(l+w 1 2 ) - 10 log(l+4w 2 ) 

(dB vs. o) 1 ,oj 2 ) (linear frequency scale) 
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Fig. 4.5a. The region of 

H(jw 1 ,ja3 2 ) = 

(l+jw 1 T 1 ) (l+jto 2 T 2 ) ( 1+ j (^2 x 2 ) 

in (log u^,log u^) plane 
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z 




z = 20 log | H | = -10 log(l+u 1 2 ) - 10 log(l+4o ) 2 ) 

2 

-10 log (I+I 6 W 2 ) 

(dB vs. log w^rlog o^) 
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lo flw. ioa. 




X-SCflLE=l . 00E+00 UNITS INCH. 
Y-SCALE=1 . OOE+OO UNITS INCH. 



Fig. 4.5c. The contour plot of the frequency response of 

. . 1 

H(]w 1 /JU) 2 ) (l+j a)l ) (l+j2u> 2 ) (l+j4a) 2 ) 

2 2 . 

z = 20 log | H | = -10 log(l+a) 1 ) - 10 log(l+4w 2 ) 

-10 log (l+16ai 2 ^) 

(log vs. log u> 2 ) 



105 




H(3u 1 /ju 2 ) (i + j Ui ) ( 1+ j 2 cj 2 ) (l+j4co 2 ) 

z = 20 log|H| = -10 log(l+u)^ 2 ) - 10 log(l + 2to 2 2 ) 

-10 log(l+16w 2 2 ) 

(dB vs. a3^,to 2 ) (linear frequency scale) 
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Fig. 4.5e. The contour plot of the frequency response of 

H(jto 1 ,ju 2 ) = (i + j aJi ) (l+j2w 2 ) ( 1+ j 4 oj 2 ) 

2 2 
z = 20 log|H| = -10 log (1+oa^ ) - 10 log(l + 4u)^ ) 

-10 log (1+16 go 2 2 ) 

(linear frequency scale) 
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Region I, when oo^ > 1/x^, 

z = apx + 20ap log x-^ dB 
Region II, when oo^ < 1/x^ 
z = 0 dB 

These separable canonic factors are identical to 
the canonic factors of a one-dimensional transfer function. 
They are just an extension to the two-dimensional case. 
Therefore, they need not be discussed further. 

B. NON-SEPARABLE FACTORS 

The following non-separable canonic factors are now 
considered . 

(1) [1 + (j<jJ 1 T 1 ) P (jw 2 T 2 ) q ] a (4 . d) 

(2) [1 + ( + ja) 2 T 2 ) P ] a (4.e) 

where p, q and a are integers. 

(1) First, consider the non-separable canonic 
factor of the form: 

H(o) 1 ,o) 2 ) = [1 + ( ju^T-L) 5 ( j« 2 x 2 ) q ] a 
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Region I, when 




z = 0 dB 



Region II, when 




2 - apx + aqy + 20ap log + 20 ay log t 2 dB 

The cornerline between the two regions is given by the 
equation 




px + qy + 20p log + 20q log t 2 = 0 

The corner line is sketched in Figure 4.6. 

The slope of the plane in region II is given by 
-a^p 2 + q 2 . Direction numbers of the normal are given 
by -ap, -aq and -1. 
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Fig . 4 . 




6. The regions of non-separable case in 
log log ^2 plane 
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Depending upon the sign of p and q and the magnitude of 
and T 2 several different orientations of the corner line 
are possible. 

Example (1) 

Discuss this canonic form in the special case with 
p = q = 1, a = +1, = 0.5, and T 2 = 1* The transfer 

function gets the following form: 

H(a) 1 ,a) 2 ) = tl + ( . 5 joi 1 ) ( joj 2 ) ] +1 

Take the logarithm of the magnitude square. 

2 4*2 

I = 1 1 ” 0 . 5 03 ^03 2 f 

z = +20 log 1 1 - . 5 to ^ 00 ^ I 

In region I, << 1 and W 2 T 2 << ^ 

z = 0 dB . 

In region II, >> 1 a nd W 2 T 2 >:> ^ 

z=x + y+20 log (0.5) 



where 

x = 20 log co 1 , y = 20 log o> 2 and z = 20 log | H (c*» 1 ,o» 2 ) 

The actual frequency response of this factor is shown in 
Fig. 4.7 including a contour plot. 



Ill 






X 






z 



T 






Fig. 4.7a. The frequency response of 

= 1 + ( j w ]_) ( j 0 • Sa^) 

2 = 20 log|H| = 20 log (1 - O.Sto^u^) 
(dB vs. log w^/log 
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'ODI *ODO\ Mfi ° 




X-SCf)LE=l . 00E + 00 UNITS INCH. 
Y-SCF!LE=1,. 05E + 00 UNITS INCH. 

Fig. 4.7b. The contour plot of the frequency response or 
H(jto 1 ,jw 2 ) = 1 + (jw L ) (j0.5ai 2 ) 
z = 20 log | H | = 20 log (1 - O.Su^u^) 
in (log co 1 ,log w 2 ) plane 
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When a = -1 we have a stability problem in this example 
The denominator of this transfer function becomes zero 
when 



1 — 0 . 5 oj ^ (a>2 — 0 



or 



This is the singularity of the transfer function in (cd-^/o^) 
plane. It can be called a singular line of the system. It 
is a curve in the plane to be compared with a pole 

in a one-dimensional transfer function. 

Example #2 : 

In this example, we discuss the same canonic form 
with different coefficients, i.e., 

p = 1 , q = -1 , a = +1 , 

t . = 0.5, T 2 = 1 

Then the transfer gets the form of 

•^ u l +i 

- U + — ) 
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Taking the logarithm of its magnitude squared 

r i 9 * ^ ^ 1 O 

| H ( to, ,ai 9 ) z = |1 + —I Z 

+ z (02 

• 5to, 

z = +20 log 1 1 + -I 

to 2 

. 5td , 

In Region I when - << 1 

“2 

z = 0 dB . 

• 5co , 

In Region II, when - >> 1 

( °2 

z = 20x - 20y + 20 log (.5) 

where x = log < 0 ^, y = log to 2 and z = 20 log |h(co^,<o 2 ) |. The 

boundary between these regions is given by the equation 

. 5to^ — to 2 

We also have a stability problem for this canonic factor when 
to 2 = 0. This is also a singularity in (io^,to 2 ) plane in the 
form of a straight line. This type of factor is unstable. 

(2) Now, consider the non-separable canonic factor 

p ct 

H((o^,oj 2 ) = [1 + (jto^to^ + j to 2 cj 2 ) ] 

We have two regions again, as follows: 

Region I, when + ^ 2 T 2 < 1 

z = 0 dB 
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Region II, when oj^x^ + a ) 2 T 2 > 1 

z = ap log (o) 1 x 1 + u) 2 x 2 ) . 



The separation line between the two regions is 
given by the equation: 



aj l T l + t0 2 x 2 = ^ 



or 



log + W 2 T 2^ = 0 . (4.6) 

The locus of points of (4.6) is shown in Figure 4.8. 
For region II, we have to think of two conditions, 
as follows: 

a) When oj^x^ > aJ 2 T 2 

z' = ctpx + ap20 log x^ 

b) When to^x-^ < oj 2 x 2 

z = aqy + ap20 log x 2 • 

Thus, region II is split into two areas with the separation 
line given by 



“l T l 



W 2 X 2 
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Fig. 4.8. Locus points of log(w^T^ + W 2 T 2^ “ 0 
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Frequency, response of this non-separable transfer function 
is shown in Fig. 4.9 with p = 1, a = -1, x = x = 1. 

C. CONCLUSIONS 

In this chapter, we discussed the properties of the 
separable and non-separable two-dimensional properties of 
canonic factors in the analog domain plane). We can 

conclude that these canonic factors can be approximated by 
planes in the log magnitude versus log and log space. 

As a design technique we can obtain any given frequency 
specification by proper combinations of planar approximation 
of canonic factors. It is obvious that only certain special 
cases and symmetries can be obtained using the factors dis- 
cussed in this thesis. For example, a wedge-shaped character- 
istic in frequency domain could not be synthesized by means 
of the factors we have considered. 

The other result of this study is the stability problem 
which arises in certain transfer functions, that is, the 
singularity curve or line. For a separable case, we don't 
have a stability problem; but in non-separable cases we may 
have stability problems. 
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z 





Fig. 4.9a. The frequency response of 

H (j“ 1 'j“ 2 ) = 1 + + jo) 2 

z = 20 log|H| = 10 log 1 1+ (w.-Ho-) | 

(dB vs. log co^,log oj ^ ) 
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0 g W 7 tOO- 



Fig . 




X -SC ALE= 1 . 00E+00 UNITS INCH. 
Y-SCALE=1 . 00E+00 UNITS INCH. 



4.9b. The contour plot of the frequency response of 

H ( joj , jw) = T 1 — — , 

2 1 + jo) 1 + 3ca 2 

2 = 20 log | H | = 10 log |l + (a^+w ) 2 | 
in (log log plane 
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z 




Fig. 4.9c. The frequency response of 

. . . , 1 

H(3 w ]_'9 w 2) 1 + + j w 2 

2 

z = 20 log | H | = 10 log |l + (u^+u^) | 

(dB vs. o) 1 ,(jo 2 ) (linear frequency scale) 



121 




v * TWO-DIMENSIONAL DIGITAL FILTER DESIGN TECHNIQUE 



A. INTRODUCTION 

The design of two and multi-dimensional recursive digital 
filters is difficult due to the absence of the Fundamental 
Theorem of Algebra in two or more dimensions [15]. That is, 
polynomials in two variables may not, in general, be factored 
into lower order polynomials. For this reason, several design 
procedures for two-dimensional recursive filters have been 
proposed in the literature. We then consider the use of the 
double bilinear z-transfer with the semi-infinite plane 
approximation in the previous chapters. These are briefly 
discussed in the following. 

Shank et al suggested three design techniques in their 
paper [11] . 

1. The use of a series of stability-preserving transforma- 
tions to produce a stable two-dimensional recursive discrete 
filters from a stable one-dimensional continuous filter. 

2. A space domain design technique in which the unit 
sample response of the filter being designed is made to 
approximate a given desired unit sample response. 

3. A stabilization technique in which the denominator 
polynomial is replaced by its double planar least-squares 
inverse. Planar least-squares can be defined as follows: 

We have some given array C. We would like to find 
an array P such that C convolved with P approximates the 
unit pulse array U. That is, 

C * P * U 
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where the symbol * denotes two-dimensional convolution. In 
general it will not be possible to make C * P exactly equal 
to U. In actuality, C * P will be equal to some other 
array, G. If we choose P such that the sum of the squares 
of the elements of U — G is minimized, P is called a planar 
least square inverse of C. Planar least square inverse of 
P is called double planar least square inverse [11] . This 
new polynomial is conjuctured to have no zeros on 
u = ^ z l ,z 2^ : i | >_, | z ^ | >_ 1} and to have a frequency 

response which closely approximates that of the original 
transfer function. 

The first technique is known as a filter design with 
rotating axes and these filters are called rotated filters. 
This design procedure consists in considering a transfer 
function in the analog frequency domain, rotating the axes, 
and applying the bilinear transformation to obtain a two- 
dimensional recursive filter. In general, the starting 
point is a separable filter and stability of the resulting 
filter is not guaranteed even if the prototype filter is 
stable. This technique was used by Costa and Venetsanopoulus 

[23] to design low-pass filters having circular symmetry, 
and also they proved that rotated filters are stable if the 
angle of rotation is between 270° and 360°. A technique 
similar in character but rotated axes in digital frequency 
domain ({z^,z^) domain) has been given by Chang and Aggarwal 

[24] . Hirona and Aggarwal [25] have accomplished the rota- 
tion through the introduction of rational powers of z^ and 
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z 2 ' an< ^ the interpolation of input and output sequences. 

They have proved that, in this technique, the resulting 
filter is stable if the prototype filter is stable. Further, 
the distortion, a consequence of the bilinear transformation, 
is absent in the resulting filter when compared with Shank's 
et al . ' s filter. 

The authors of [11] were unable to prove that the new 
denominator polynomial produced by the third technique , 
which is an extension of the one-dimensional case, was 
stable in the two-dimensional case. Nonetheless, several 
hundred filters have been designed using this technique 
without a counter example appearing [18]; furthermore. Jury, 
Kolavennu and Anderson [26] have shown it to be valid for 
certain low order cases. However, Genin and Kamp [27] have 
recently shown that counter-examples do, in fact, exist. 

Spectral factorization is another technique which can be 
used to produce a stable filter from an unstable one. By 
spectral factorization we mean a complex map that carries a 
stable rational transfer function into another stable transfer 
function exhibiting a different frequency response, at the 
same time maintaining some desirable characteristic [29]. 
Pendergrass, Mitra and Jury [28] used stability-preserving 
spectral transformations to modify the frequency response of 
quarter-plane filters. 

Mitra and Chakrabarti [29] extended this work. Read 
and Treitel [30] suggested a two-dimensional Hilbert trans- 
form to perform the spectral factorization, while Pister 
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4 










[31] and Ekstrom and Woods [32] made use of the two- 
dimensional cepstrum. Ekstrom and Woods also suggested ways 
that the factors could be truncated and smoothed to preserve 
stability. They also showed that the usual quarter-plane 
filter is not the most general two-dimensional recursive 
filter. It is instead possible to have a recursive filter 
which in the limit has a unit sample response with support 
on an asymmetric half-plane, and also they developed 
stability conditions for this new class of filters. 



B. TWO-DIMENSIONAL DIGITAL FILTER DESIGN PROCEDURE BY 
DOUBLE Z-TRANSFORM 

1. Double Z-transform 

In the one-dimensional recursive digital filter 
design, the bilinear z-transform is often applied to an 
analog filter transfer function. We proposed that this one- 
dimensional design approach can be extended to the two- 
dimensional case. Before going into the design procedure, 
two-dimensional recursive filter design via the double 
bilinear z-transform, it is better to prove that we can 
use double z-transform for two-dimensional systems. 

Each of the canonic factors, proposed in Chapter 
IV, either separable or non-separable, can be converted 
into discrete equations by the double z-transform, as 
follows : 



S 



1 





(5.1a) 
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(5.1b) 



2 ^ Z 2 
" 2 = ^ ( rrz i) 



The single z-transform is equivalent to trapezoidal 
integration which can be shown as follows: 



X (s) = S Y (s ) 



(5.2) 



x(t) = 



dy(t) 

dt 



or 



Y 



n 



Yn- 



n-1 



n 

/ x(t) dt 

n-1 



T, 

2 (x n 



+ x n-l } 



(5.4) 



where T is the sampling period. 

We write (5.4) in the z transform domain as 

(1 - z” 1 ) Y (z ) = |(l + z _1 )X(z) 



Hence 



X ( z ) 





+ z 



Y(z) 



(5.5) 
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Comparing (5.5) with (5.2), it is seen that 



S 



1(1_Z 
T 1 + 




(5.6) 



Equation (5.6) is known as a trapezoidal bilinear z-transform. 

We can drive the double z-transform in a similar 
way, as follows: 



W(x,y) 




dxdy 



Z(x,y) 



(5.7a) 



or 



W(s^,s 2 ) = s i s 2 z ( s i' s 2^ 



(5.7b) 



Then 



/ / d 2 Z(x,y) = / / W(x,y) dx, dy (5.8) 

by a discrete approximation, the left-hand side of (5.8) 
becomes (see Fig. 5.1) 

/ / d 2 Z(x,y) = [Z(x n ,y n - Zix n 'Y n -i)l 

~ IZ^n-l'V ' Z(x n-l' y n-l )] 

= (1 - z" 1 - z' 1 + z^z” 1 ) Z (z-^z^ (5,9) 
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and the right-hand side of (5.8) is approximated as follows: 



//w(x,y) dxdy = 



^[W(x n ,y n ) + wu^y^) 



+ M(x n-l' y n> + W(x n-l' y n-l )] 



AxAy ri , -1 , -1 , -1 -l 1rT/ , 

= — 4 + z i + z 2 + Z 1 z 2 ] w ( z ]_' z 2^ 



By equation (5.10) and (5.9) we have 



(5.10) 



W(z 1 /Z 2 ) 1 - z” 1 - z" 1 + z^z' 1 4 

Z(z l' z 2 } = 1 + z" 1 + z" 1 + z^z" 1 ‘ **** 



(5.11) 



or 



i _1 i _1 

2 1 ~ z i 2 ^ ” z 2 



1 + z. 



1 + z 



(5.12) 



By comparing (5.1a) and (5.2b) it is seen that 



2 1 z 

aI ( — ^ 

1 + z^ 

2 1 ' Z 2 1 

^rrz 1 ' 



Q . E . D . 



Thus, this is called a double z-transform. The 
double z-transform seems to be a reasonable substitution 
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yO) 


d 


X (O 




d t 






(a) One-dimensional case 



Z( x »y) 


d 2 


W(x,y) 




d x dy 






Fig. 5.1. Two-dimensional extension of the 
Bilinear Transformation 
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for each canonic factor. However, each factor must be 
checked for stability, which is discussed in detail later in 
this chapter. If each transformed factor is BIBO stable the 
overall transfer function will be BIBO stable. But sometimes the 
double z-transform may in certain cases lead to unstable 
solutions [33] . Finally, it is noted that the double z-transform 
results in a quarter-plane symmetry which is a severe limitation. 

2 . Design Procedure 

The two-dimensional recursive digital filter design 
investigated in this thesis is the double bilinear z-transform. 
Proceed as follows : 

1. Determine a two-dimensional analog transfer 
function. 



H (s-j^, s 2 ) 



A(s 1 ,s 2 ) 
B (s~,s 2 ) 



(5.12) 



where A(s lf s 2 ) and B(s 1 ,s 2 ) are mutually prime two-dimensional 
polynomials and B(s^,s 2 ) satisfies 



B(s 1 ,s 2 ) 0 

(s^' s 2) £ ^ = ( ( s i' s 2 ) * 



(5.13) 

Re (s 1 ) >. 0,Re(s 2 ) ^ 0} 



in such a way that, the frequency response of G(s 1 ,s 2 ) 
meets the given specifications (for example as discussed 
in Chapters III and IV) . 
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2. Apply the double bilinear z-transform 



S 



1 



S 



2 







(5.14a) 



(5.14b) 



where A is the sampling period, which is assumed the same in both 
directions. This yields the two-dimensional digital transfer 
function 



2 ^ Z 1 

H 'i<— ^ t> 

1 + Z~ 





1 + z' 1 



) ) 



T(z r z 2 ) = 



C (Z]_ ' z 2 ) 
D ( z ]_ ' z 2 ) 



(5.15) 



where C(z^,z 2 ) and D(z^,z 2 ) are mutually prime two-dimensional 
polynomials in z^ and z 2 « (5.15) will have the desired 

frequency response with proper sampling time without 
aliasing effect. Since 



2 1 - Z i 1 

s i - 



maps the closed unit disc onto the closed right half 
plane, one would expect (5.15) to satisfy the following 
sufficient stability conditions [11], [20]. 



D(z 1 ,z 2 ) ? 0 
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{ ( 2 1 ,z 2 ) : IzJ > 1, |z 2 | > 1} 



(5.16) 



\/ ( 2 x / z 2 ) e u 2 = 



Unfortunately, however, the conditions (5.13) and (5.16) 
as we will show later in this chapter, are not always 
satisfied . 

In summary , we have two conditions to be satisfied 
for stability. 

I. B (s 1 ,s 2 ) ^ 0 

\/( s 1 ' s 2 ) e r = {(s 1 ,S 2 )s R e s 1 _> 0,R e s 2 >_ 0} (5.17) 

II. D(z 1 ,2 2 ) fi 0 

A 

Y(z 1 , 2 2 ) £ u = { (z 1 ,z 2 ) : | z ± | > 1, | z 2 | >_ 1} 

(5.18) 

where V stands for "for each", and u 2 is the closed unit bidisc. 

C. DIGITAL CANONIC FACTORS AND STABILITY CONSIDERATIONS 
1 . Introduction 

Certain properties of two variable polynomials and 
rational functions which will be needed later will now 
be discussed. It is well known that a two variable polynomial 
is not, in general, factorable into first-order polynomials; 
but, a two variable polynomial can be factored into irreduci- 
ble factors which are themselves two variable polynomials 
but which cannot be further factored. Of course a given 
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polynomial may itself be irreducible. These irreducible 
polynomials are unique up to multiplicative constants [33] . 
Two polynomials which have no irreducible factors in common 
are said to be mutually prime. Consider two variable 
rational functions, in general. 



G(z 1 ,z 2 ) 



P (z 1 ,z 2 ) 
Q(z 1 ,z 2 ) 



where P(z^,z 2 ) and Q(z^,z 2 ) are mutually prime. Using 
terminology of [34], a point (z^,z 2 ) making Q(z£,zJ,) = 0 
but P(z^,z 2 ) ^ 0 is called a pole or a nonessential singu- 
larity of the first kind (such a point is analogous to a 
pole in the one variable case). A point (z£,z 2 ) making 
P(z£,z 2 ) = Q(z£,z 2 ) = 0 is called a nonessential singu- 
larity of the second kind (such a point has no analogue in 
the one variable case) . Clearly, if (z£,z 2 ) is a pole, 
G(z^,z 2 ) = », but if (zj,zj) is a nonessential singularity 
of the second kind, then G(z£,z 2 ) is undefined. 

2. Separable Factors 

The most common canonic separable factor in Chapter 

IV is 



H (o)^ , cu 2 ) 



(1 + joj^T^) (1 + j(i) 9 T 9 ) 



(5.19a) 



2 2 ; 



assuming = t 2 = !• Substituting 



jo-, — S-, , 
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j“2 “ S 2 

We get 



H(S 1 ,S 2 ) " (1 + sjd + sJ (5.19b) 

Substituting equations (5.1a and 5.1b) into (5.19b), 
perforins the double z-trans formation. Thus, 

S (1 + z" 1 ) (1 + z" 1 ) 

T(z x ,z 2 ) = Kr £— (5.20) 

(a + z 1 ) (a + z^) 



where 



B 



(A ~ 2) 



2 



a 



a + 2 
a — 2 



A = sampling period 

This separable canonic form obviously satisfies the condition 
of (5.17). We should be careful about the second condition 
of (5.18) . 

When z i = z 2 ~ we have a nonessential singularity 

of the second kind if a = 1, but a = 1 never occurs. 
Therefore, there is no nonessential singularity of the second 
kind. Thus, this separable case is stable by satisfying 
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two essential necessities (5.17) and (5.18) . The frequency 
response of this separable case is plotted in Fig. 5.2 , with 
sampling time A = 0.08 which fits that of the analog frequency 
response . 

The double z-transform for equation (5.19a) is 



T * Z 1' Z 2* 



(A- 2 t^) (A^t^) A +2 



(1 + z" 1 ) (1 + z" 1 ) 



A +2 



, T i , _-i x , 

( A- 2 t 1 Z 1 } A -2 t 2 



t 2 , - 1 . 
+ z 2 ) 



(5.21) 



and this general case for the two-dimensional separable 
factor is stable by satisfying the stability necessities 



of (5.17) and (5.18) 



The second analog canonic factor in Chapter IV is 



H(j Ul /ju 2 ) (1 + ja) 1 T 1 ) (1 + jcu 2 x 2 ) (1 + jo> 2 T 3 ) * 



When we perform the double bilinear transform for 
this, we get 



T(z l' z 2 ) 



3 -1 -12 

A J (1 + z 1 X ) (1 + z 2 V 



A +2 x -i A +2 

(A- 2 x 1 ) (A- 2 t 2 ) (A- 2 t 3 ) [ a _ 2t +z i 1 ^a- 2 t 



T 2 . -1. r 
+ z 2 ][ 



A+2t, -I 

-- ,- - + zl x ] 

A-2t 2 2 



(5.22) 



The frequency response of (5.22) for =1, t 2 = 2, = 4 

is shown in Fig. 5.3 which almost fits the analog case. 
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z 











Fig. 5.2a. The frequency response of 

3 (1+z, -1 ) (1+z " X ) 

T( z z ) = ^ ^- T - 

(a+z x x ) (a+z 2 x ) 

2 

where 6 = — j , a = , A = 0.0 8 

(A-2) Z a Z 

(dB vs. log log cd^) 
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z 





Fig. 5.2b. The frequency response of 



T (z^ , z^) 



6(l+z 1 1 ) ( 1+z 2 1 ) 
(a+z.^ 1 ) (a+z 2 -1 ) 



6 = 



(A-2) 



2 ' 



a = , A = 0.08 



(dB vs. (linear frequency scale) 
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Fig. 5.2c. 



The contour plot of the frequency response of 
B(l+z 1 -1 ) (l+z 2 -1 ) 



T(z 1 ,z 2 ) = 



e = 



2 ' 



(a+z 1 “ 1 ) (a+z 2 1 ) 

a = |±| , A = 0.08 



where 



(A-2) 

(linear frequency scale) 
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We can summarize the m xn analog separable factor as 



H ( j <D^ , j ( 1)2 ) 




(5.23) 



After performing double z-transform, (5.23) becomes 



T(2 1 ,z 2 ) 



M 
[ n 
i=l 




] [ n 

j=i 



N 




1 5 .24) 

Equations (5.23) and (5.24) satisfy the conditions 
of (5.17) and (5.18). As a result, a separable filter 
results in separable filters in digital domain and they are 
stable . 

3 . Non-separable Factors 

The first non-separable factor in Chapter IV is 
4 . d with a = - 1 , p = q = x 1 = t 2 = 1 as given below 



H ( j o)^ , j ^2 ) 



1 



(5.25a) 



1 4" J (i] ^ J cp^ 



or 



H (s 1 ,s 2 ) 



1 



(5.25b) 
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X 



z 

_L 



T 




where A = 0.08 
(dB vs. log , log u^) 



r 
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z 





T( Zl ,z 2 ) = 



A 3 (l+z 1 ~ 1 ) (l+z 2 _1 ) 



/A T\ /A A \ f A 0\ / A + 2 —1. , A+4 -1 x , A+8 — 1 x 

(A-2) (A 4) (A 8) (^_ 2 +z x ) +z- L ) +z 2 ^ 



where A = 0.08 



(dB vs. u)^,co 2 ) 
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' 0 . 



0 - 



5 . 



! 0 . 




Fig. 5.3c. The contour plot of the frequency response of 



T(z r z 2 } = 



A 3 (l+z 1 1 ) (l+z 2 _1 ) 



(A-2) ( A — 4 ) (A-8) (f±| tz^ 1 ) (f±i tz^ 1 ) (|±| tz^ 1 ) 



where A = 0.08 
(linear frequency scale) 
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and by the double z-transform, we get 



T(z lf z 2 ) = 



A 2 (l +z 1 1 ) (1 + z" 1 ) 



(A 2 + 4) (a 2 -4) (z ” 1 + z" 1 ) + (a 2 + 4)z^ 1 z 2 1 



(5.26) 



In the analog domain (5.25b) does not satisfy the 
condition (5.17). We have a singularity at 



u l u 2 



1 



For the second condition, we have a nonessential 
singularity of the second kind at z^ = -1, z n = -1. This 
does not tell us whether it is stable or not. When we 
check the frequency response, it is seen that this factor 
is unstable. 

The second non-separable canonic factor can be 
derived from (4.d) with a = q = -1 and = T 2 = 1 as 
follows . 



H ( j , j 0 ) 2 ) 



1 + 



3 oj- 
j W » 



(5.27a) 



or 



H (s 1 , s 2 ) 




S 1 + s 2 



(5.27b) 



14 3 



Then, (5.27b) becomes after performing the double z-transform, 
as follows: 



T ( z ]_ ' z 2 ) 



(1 + z 1 1 ) (1 - z 2 X ) 
2(1 - z" 1 ^ 1 ) 



(5.28) 



This canonic factor does not satisfy either 
stability condition. Hence it is unstable. 

The other series of canonic factors can be derived 
from (4 . e) . 

For example, with a = -1, t 1 =t 2 =p= 1, we 

get 



H ( j w 



1/ jco 2 ) 



1 + jio-, + jw- 



(5.29a) 



or 



H ( s , s 2 ) 



1 




or in the digital domain 



(5.29b) 



A (1 + z, 1 ) (1 + z- 1 ) 

T(Z.,Z ) = ry ^ (5.30) 

(A+4) +A(z 1 i + z 2 ) + (A-4)z 1 i z 2 i 

This canonic factor satisfies the first condition but 
it has a nonessential singularity of the second kind, at 
z^ = z 2 = -1 and so we should check the frequency 
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response. We see that there is a singularity at 

^ C °1 ,C0 2^ = ( ^ / tt ) • Thus this non-separable canonic factor 
is unstable. 

4 . Conclusions 

All of the canonic factors discussed in this 
chapter are collected in Table I. By comparing the results 
shown in this table, we can say that separable stable canonic 
factors result in stable separable filters. As seen from 
Table I, the separable transfer functions 1, 2 and 3 which 
start with stable separable ones in the two-dimensional 
continuous frequency domain end up as stable separable trans- 
fer functions in the two-dimensional discrete domain. This 
is also true for the mxn separable canonic factor 4. 

Non-separable canonic factors have very different 
characteristics. The unstable non-separable case ends up 
with a non-separable case in the discrete domain, which is 
also unstable. Although we could not prove this result 
generally, we also could not find any counter example to 
exist. The transfer functions of 5, and 6 start in the 
unstable case, and end up with the unstable case in the 
discrete domain. The transfer functions of 7 starts with 
the stable case in the continuous frequency domain, ends up 
with the unstable case in the discrete domain. Thus we 
have to be very careful with the double z-transform since it 
does not yield stable results every time. It may yield 
unstable results even if it starts from the stable case. 
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VI. CONCLUSIONS 



A technique for the use of semi-infinite planes to 
approximate the log magnitude versus log frequency charac- 
teristics of two-dimensional filters is presented. This 
technique is applied to quarter plane filters and works 
well for separable transfer functions. Other non- 
separable canonic (basic) transfer functions are also 
studied in terms of their planar approximation. 

The following separable and nonseparable canonic factors 
have been studied. 

(1) (1 + jto 1 x 1 ) p (l + jto 2 T 2 ) g 

(2) (1 + to^T^) p (l + j w 2 T 2 ) ^ ( 1 + 

(3) [1 + (jw 1 T 1 ) p ] a 

(4) [1 + (jt0 1 T 1 ) P (l + j(0 2 T 2 ) g ] a 

(5) [1 + (jw 1 T 1 + jco 2 x 2 ) p ] a 

The properties and stability of these canonic factors 
are investigated first in the continuous variable domain 
and then after the double bilinear z-transform is applied 
to extend them to two-dimensional recursive digital filters. 
These are all summarized in Table I in Chapter V. An 
important result is that although the double bilinear 
z-transform works well for separable factors, it does not 
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always lead to a stable digital realization starting with a 
stable analog transfer function. 

The following are suggested for further research: 

1. Other types of canonic factors need to be considered. 
For example, the general form: 

H(j03 1 ,j0J 2 ) = [1 + (jc0 1 T 1 ) P +(jw 2 T 2 ) q +(jc0 1 T 3 ) r (j(0 2 X 4 ) t ] a 

2. The work should be extended to non-causal (half- 
plane) filters. 

3. Transfer functions which have circular symmetry in 
the discrete domain need to be investigated. 

4. Another approach is to consider half-plane 
recursive algorithms, by working backwards from the discrete 
domain (z^,z 2 ) to continuous domain (s^,s 2 ). 

5. Stability of canonic factors needs to be investigated 
in detail. 

6. The use of predistortion should also be considered. 
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